An introduction to `patter`

Edward Lavender

edward.lavender@eawag.ch

Eawag, Swiss Federal Institute for Aquatic Science and Technology

Introduction
Methodology & outputs
Set up
Workflow
- Overview
- Boilerplate setup
- An example with real data
- Simulation-based workflow
Ancillary functions
Advanced topics
Development and help

Introduction

https://github.com/edwardlavender/patter

Introduction

https://github.com/edwardlavender/patter-workshops

Introduction

patter is a R package in the movement-ecology ecosystem
patter was motivated by our research in passive acoustic telemetry systems
patter fits state-space models to animal-tracking data, using particle filtering and smoothing algorithms

Team

Institutions

Workshop

The goal today is to introduce patter
We will focus on how to use the package
There will be:
- Interactive presentation & questions
- Optional practical activity (using example datasets or your own)

Questionnaire

Who has managed to install patter and connect it to Julia?
Experience with animal-tracking?
Experience with R and Julia?
Experience with statistical modelling, including SSMs?
Experience with Bayesian statistics?
Experience with particle algorithms?

Context

There are many packages for animal tracking (Joo et al., 2020)
Numerous packages are available for passive acoustic telemetry (VTrack, glatos, RSP, etc.)
patter is a unique addition that fits SSMs to animal-tracking data using particle and smoothing algorithms
patter can incorporate a wide range of data types

Important

For passive acoustic telemetry, patter is the only package that reconstructs individual movements and patterns of space use within a coherent probabilistic framework.

Context

patter uses a dedicated high-performance backend written in Julia

https://github.com/edwardlavender/Patter.jl

Evolution

patter evolved from flapper (Lavender et al., 2023)

Evolution

flapper is a more general-purpose package for passive acoustic telemetry that supports:
- data acquisition, simulation, assembly and processing
- spatial operations and distance metrics (e.g., fast C++ LCP algorithms)
- data exploration (detection statistics, movement metrics)
- ‘flapper’ algorithms

Note

Let us know if there are routines in flapper that you’d like to see in patter.

Evolution

patter focuses specifically on inference for SSMs using particle algorithms

https://github.com/edwardlavender/patter

Methodology

For the methodology, see:

Lavender, E. et al. (2025a). Particle algorithms for animal movement modelling in receiver arrays. Meth. Ecol. Evol. https://doi.org/10.1111/2041-210X.70028

For the package, see the package paper:

Lavender, E. et al. (2025b). patter: particle algorithms for animal tracking in R and Julia. Meth. Ecol. Evol. https://doi.org/10.1111/2041-210X.70029

For a case study, see:

Lavender, E. et al. (2025c). Animal tracking with particle algorithms for conservation. BioRxiv. https://doi.org/10.1101/2025.02.13.638042

Methodology

For a algorithmic explanation of the methodology for mathematical and non-mathematical readers:

Lavender, E. et al. (2025a). Supporting Information.

For an intuitive summary, see:

vignette("a-methodology", package = "patter")

Methodology

We formulate a Bayesian state-space model for the joint probability distribution $f(\mathbf{s}_{1:T} \mid \mathbf{y}_{1:T})$ of a tagged individual’s possible trajectories from the start to the end of the time series
The joint probability distribution is expressed in terms of a prior (the movement process) and a likelihood (the observation process)

\[f(\mathbf{s}_{1:T} \mid \mathbf{y}_{1:T}) \propto f(\mathbf{s}_{1:T}) f(\mathbf{y}_{1:T} \mid \mathbf{s}_{1:T})\]

Methodology

Inference for the latent states is performed using particle filtering and smoothing
The particle filter approximates the partial marginal distribution $f(\mathbf{s}_t \mid \mathbf{y}_{1:t})$ using a set of weighted samples called ‘particles’
The particle smoother approximates the full marginal distribution $f(\mathbf{s}_t \mid \mathbf{y}_{1:T})$

Methodology

Particles are like grains of sand

Methodology

A particle filtering algorithm

Methodology

Particle smoothers reweight particles from the filter to approximate $f(\mathbf{s}_{t} \mid \mathbf{y}_{1:T})$
The two-filter smoother consists of a re-weighting of particles from a forward filter run and a backward filter run based on the probability density of movements between pairs of locations

Methodology

The two filter smoothing algorithm

Methodology

Smoothing is like pruning

Outputs

Why should I use `patter`?

patter gives us a probabilistic representation of the possible locations of an individual though time that incorporates our ecological knowledge and all the data that we collect

Why should I use `patter`?

Refined maps of space use:
- with a proper biological and statistical interpretation
- which outperform those from heuristic models across the board
Refined estimates of residency:
- At time-scales of interest (e.g., daily, weekly, yearly residency)
Refined downstream analyses (e.g., habitat preferences) that can incorporate uncertainty in the individual’s location through time

Why should I use `patter`?

patter produces refined maps of space use (Lavender et al., 2025a)

Why should I use `patter`?

patter consistently outperforms heuristic methods (Lavender et al., 2025a)

Why should I use `patter`?

patter produces refined estimates of residency (Lavender et al., 2025b)

Set up

Operating systems

patter v.1.0.0 supports Windows and MacOS
patter v.2.0.0 supports Windows, MacOS and Linux

Note

There are some special considerations for running patter on linux that we won’t discuss today. Please get in touch if you are running patter on HPC linux environments.

Installation

To install patter, follow the instructions on https://github.com/edwardlavender/patter
- Install R and Julia
- On Windows, install RTools
- On MacOS, if required, install Xcode (and configure .Makevars)
- Install patter (and let it handle Julia dependencies)

Important

A recent Julia version (≥ 1.10 [LTS]) is recommended. Download links are here or you can use juliaup. Follow the README instructions for further guidance.

Project set up

(optional) Set up an RStudio project
(optional) Initiate local-dependency management with renv (use renv::init())
(optional) Build project directories (data/, data-raw/, R/, Julia/, …)
Install, load & attach patter (as explained previously)
(optional) Set Julia options (next slide)
Connect to Julia (via julia_connect())

Note

renv is an R package that implements local-dependency management. In Julia, we get this ‘for free’. In patter, set JULIA_PROJ to use a local package environment for the Julia dependencies (this is generally recommended).

Project set up: `JULIA` options

JULIA_HOME is the location of the Julia installation
- Usually, JULIA_HOME is found automatically, but it may be required if Julia is installed in a non-standard location
- E.g., on an Intel Mac JULIA_HOME may be something like /Applications/Julia-1.10.app/Contents/Resources/julia/bin/
(optional) JULIA_PROJ is the location of the Julia project (./Julia/)
(optional) JULIA_NUM_THREADS is the number of threads
(optional) JULIA_PATTER_SOURCE is the source of the Patter.jl package (for advanced use)

Important

On Windows, JULIA_NUM_THREADS must be set system-wide as an environment variable. Follow the guidance here.

Project set up: `JULIA` options

It is a good idea to set JULIA options globally via .Rprofile (or .Renviron):

# Set JULIA_HOME
# * This is usually not required

# (optional) Set JULIA_PROJ to ./Julia/
if (!requireNamespace("here", quietly = TRUE)) {
  install.packages("here")
}
Sys.setenv(JULIA_PROJ = here::here("Julia"))

# (optional) Set JULIA_NUM_THREADS
# * This works on MacOS and Linux 
Sys.setenv(JULIA_NUM_THREADS = parallel::detectCores())

# (optional) Set additional options 
Sys.setenv(OMP_NUM_THREADS = parallel::detectCores())

Note

See ?julia_connect() for guidance.

Connect to `Julia`

After library(patter), run julia_connect() to connect R to Julia. This:
1. Sets up JuliaCall (which handles the R–Julia coupling)
2. Validates the Julia installation (@ JULIA_HOME)
3. (Optional) Activates a local Julia project (JULIA_PROJ)
4. Installs and pre-compiles Julia dependencies (e.g., Patter.jl)
In general, you should run julia_connect() once per R session.

Note

The first call to julia_connect() may take several minutes (or longer) as Julia packages are installed and pre-compiled. (On old laptops, we have observed times up to 30 minutes!) Subsequent calls are faster.

Connect to `Julia`

The syntax for julia_connect() is:

julia_connect(JULIA_HOME, JULIA_PROJ, JULIA_NUM_THREADS, ...)

If you’ve set environmental variables (recommended), just use:

julia_connect()

Workflow

The patter workflow

Workflow

We will work though:
- Overview
- Boilerplate setup
- An example workflow with real data
- An example simulation-based workflow

Note

For a summary of the workflow, see vignette("b-workflow-outline", package = "patter"). See the README for an example workflow. patter functions are well documented. Please read them carefully and follow the links to the Patter.jl documentation where directed. Help us to improve the documentation by submitting GitHub issues.

Workflow

patter is easy:
1. Define the state type and the movement model
2. Define the observation model(s)
3. Simulate observations or assemble real-world datasets
4. Run the forward filter and the backward information filter
5. Run the particle smoother
6. Analyse the results!

Workflow

The challenge is to formulate biologically sensible movement and observation models
The code is just boilerplate: copy & paste and edit as required!
Please reach out on GitHub with queries

Overview

Boilerplate setup (1/7)

# Load & attach package 
library(patter)
op <- options(patter.verbose = FALSE)

# Ancillary packages
library(data.table)
library(dtplyr)
library(dplyr, warn.conflicts = FALSE)
library(glue)
library(JuliaCall)
library(spatstat.explore)
library(truncdist)

Boilerplate setup (1/7)

# Load datasets
head(dat_detections, 2)

   individual_id           timestamp receiver_id
           <int>              <POSc>       <int>
1:            25 2016-03-17 01:50:00          26
2:            25 2016-03-17 01:52:00          26

head(dat_moorings, 2)

   receiver_id receiver_start receiver_end receiver_x receiver_y receiver_alpha
         <int>         <POSc>       <POSc>      <num>      <num>          <num>
1:           3     2016-03-03   2016-11-02   706442.1    6254007              4
2:           4     2016-03-03   2017-03-08   709742.1    6267707              4
   receiver_beta receiver_gamma
           <num>          <num>
1:         -0.01            750
2:         -0.01            750

head(dat_archival, 2)

   individual_id           timestamp  depth
           <num>              <POSc>  <num>
1:            25 2016-03-16 00:00:00 159.15
2:            25 2016-03-16 00:02:00 161.01

Boilerplate setup (1/7)

# Connect to Julia 
julia_connect()
set_seed()

Boilerplate setup (1/7)

# Set map 
map <- dat_gebco()
set_map(map)

# Study timeline 
timeline <- seq(as.POSIXct("2016-03-17 01:50:00", tz = "UTC"),
                as.POSIXct("2016-03-18 01:48:00", tz = "UTC"), 
                by = "2 mins")

Set state (2/7)

state <- "StateXY"

Set movement model (3/7)

mobility   <- 750.0
model_move <- model_move_xy(.mobility = mobility, 
                            .dbn_length = glue("truncated(Gamma(1.0, 250.0), upper = {mobility})"),
                            .dbn_heading = "Uniform(-pi, pi)")
plot(model_move)

Set observation model(s) (4/7)

# Acoustic observations 
det <- dat_detections[individual_id == 25L, ]
acc <- assemble_acoustics(.timeline   = timeline,
                          .detections = det,
                          .moorings   = dat_moorings, 
                          .services   = NULL)
head(acc, 3)

             timestamp sensor_id   obs receiver_x receiver_y receiver_alpha
                <POSc>     <int> <int>      <num>      <num>          <num>
1: 2016-03-17 01:50:00         3     0   706442.1    6254007              4
2: 2016-03-17 01:50:00         4     0   709742.1    6267707              4
3: 2016-03-17 01:50:00         7     0   708742.1    6269107              4
   receiver_beta receiver_gamma
           <num>          <num>
1:         -0.01            750
2:         -0.01            750
3:         -0.01            750

Set observation model(s) (4/7)

# Archival observations 
arc <- dat_archival[individual_id == 25L, ]
arc <- assemble_custom(.timeline = timeline, 
                       .dataset = 
                           arc |> 
                             mutate(sensor_id = 1L) |> 
                             select("sensor_id", "timestamp", obs = "depth") |>
                             mutate(depth_sigma = 50, depth_deep_eps = 20) |>
                           as.data.table())
head(arc, 3)

   sensor_id           timestamp   obs depth_sigma depth_deep_eps
       <int>              <POSc> <num>       <num>          <num>
1:         1 2016-03-17 01:50:00 73.78          50             20
2:         1 2016-03-17 01:52:00 73.32          50             20
3:         1 2016-03-17 01:54:00 73.32          50             20

Set observation model(s) (4/7)

# All observations
yobs <- list(ModelObsAcousticLogisTrunc = acc, 
             ModelObsDepthNormalTruncSeabed = arc)

Run particle filter (5/7)

# List filter arguments
args <- list(.timeline    = timeline,
             .state       = state,
             .yobs        = yobs,
             .model_move  = model_move,
             .n_particle  = 1e5L)

# Forward run
fwd <- do.call(pf_filter, args, quote = TRUE)

# Backward run
args$.direction <- "backward"
bwd <- do.call(pf_filter, args, quote = TRUE)

Run particle smoother (6/7)

smo <- pf_smoother_two_filter()

Analyse outputs (7/7)

head(smo$states)

   path_id timestep           timestamp map_value        x       y
     <int>    <int>              <POSc>     <num>    <num>   <num>
1:       1        1 2016-03-17 01:50:00  84.47292 709203.3 6253549
2:       1        2 2016-03-17 01:52:00  66.52094 709153.5 6253109
3:       1        3 2016-03-17 01:54:00 108.30474 709063.3 6253516
4:       1        4 2016-03-17 01:56:00  96.38661 709127.0 6253461
5:       1        5 2016-03-17 01:58:00  66.52094 709181.6 6253066
6:       1        6 2016-03-17 02:00:00  84.47292 709198.5 6253517

Analyse outputs (7/7)

head(smo$diagnostics)

   timestep           timestamp       ess maxlp
      <int>              <POSc>     <num> <num>
1:        1 2016-03-17 01:50:00 1000.0000   NaN
2:        2 2016-03-17 01:52:00  838.4857   NaN
3:        3 2016-03-17 01:54:00  827.4884   NaN
4:        4 2016-03-17 01:56:00  597.2746   NaN
5:        5 2016-03-17 01:58:00  851.9332   NaN
6:        6 2016-03-17 02:00:00  821.2493   NaN

Analyse outputs (7/7)

smo$callstats

             timestamp              routine n_particle n_iter loglik
                <POSc>               <char>      <int>  <int>  <num>
1: 2025-06-01 13:48:58 smoother: two-filter       1000     NA    NaN
   convergence     time
        <lgcl>    <num>
1:        TRUE 6.047296

Analyse outputs (7/7)

map_dens(.map = map, .coord = smo$states, .sigma =  bw.diggle)$ud

class       : SpatRaster 
dimensions  : 264, 190, 1  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 695492.1, 714492.1, 6246657, 6273057  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 29N (EPSG:32629) 
source(s)   : memory
name        :         lyr.1 
min value   : -3.217162e-19 
max value   :  3.654440e-03

Boilerplate setup

Load and attach patter and supporting packages:

# Load & attach patter
library(patter)

# Ancillary packages
library(data.table)
library(dtplyr)
library(dplyr, warn.conflicts = FALSE)
library(glue)
library(JuliaCall)
library(spatstat.explore)
library(truncdist)

Use patter.verbose to control user output messages:

# (patter.verbose = TRUE by default)
op <- options(patter.verbose = FALSE)

See ?patter-progress for additional options

Boilerplate setup

Connect to Julia:

# Connect to Julia
# (Specify `JULIA_PROJ` and `JULIA_NUM_THREADS` in .Rprofile)
julia_connect()

Set the seed in R and Julia via set_seed():

set_seed(2024)

An example workflow with real data

Study system

This first part of this presentation walks though how to reconstruct movements and patterns of space use using example data collected from flapper (Dipturus intermedius) in Scotland.

The Loch Sunart to the Sound of Jura Marine Protected Area

Study system

Read observations
Define map of study system
Define timeline of interest
Assemble observations

Study system

Observations include detections at receivers:

head(dat_detections)

   individual_id           timestamp receiver_id
           <int>              <POSc>       <int>
1:            25 2016-03-17 01:50:00          26
2:            25 2016-03-17 01:52:00          26
3:            25 2016-03-17 01:54:00          26
4:            25 2016-03-17 01:58:00          26
5:            25 2016-03-17 02:00:00          26
6:            25 2016-03-17 02:04:00          26

Note

patter uses the words detections and acoustics to distinguish detection data (which comprise detections alone) and acoustic data (which comprise detections and non-detections, which are implicitly known for each time step).

Study system

This is the receiver metadata:

head(dat_moorings)

   receiver_id receiver_start receiver_end receiver_x receiver_y receiver_alpha
         <int>         <POSc>       <POSc>      <num>      <num>          <num>
1:           3     2016-03-03   2016-11-02   706442.1    6254007              4
2:           4     2016-03-03   2017-03-08   709742.1    6267707              4
3:           7     2016-03-03   2016-11-26   708742.1    6269107              4
4:           9     2016-03-03   2016-10-13   706042.1    6254307              4
5:          11     2016-03-03   2016-11-18   707542.1    6267707              4
6:          12     2016-03-03   2017-03-08   710042.1    6267307              4
   receiver_beta receiver_gamma
           <num>          <num>
1:         -0.01            750
2:         -0.01            750
3:         -0.01            750
4:         -0.01            750
5:         -0.01            750
6:         -0.01            750

Tip

To incorporate passive acoustic telemetry datasets in patter, follow the format for dat_detections and dat_moorings. There is a pat_setup_data() function that you can use to check these datasets. The receiver_alpha, receiver_beta and receiver_gamma columns optionally define receiver-specific parameters of a truncated, distance-decaying logistic detection probability model but (as we will see later) we can use any observational model of our choosing.

Study system

We also have archival (depth) observations, collected at a resolution of two-minutes:

head(dat_archival)

   individual_id           timestamp  depth
           <num>              <POSc>  <num>
1:            25 2016-03-16 00:00:00 159.15
2:            25 2016-03-16 00:02:00 161.01
3:            25 2016-03-16 00:04:00 162.86
4:            25 2016-03-16 00:06:00 164.48
5:            25 2016-03-16 00:08:00 158.46
6:            25 2016-03-16 00:10:00 155.45

Study system

Our goal in the first part of this presentation is to reconstruct movements and patterns of space use for an example individual given these data.

# Define detections
det <-
  dat_detections |>
  filter(individual_id == 25L) |>
  mutate(individual_id = NULL) |>
  as.data.table()

# Define archival (depth) observations
arc <-
  dat_archival |>
  filter(individual_id == 25L) |>
  mutate(individual_id = NULL) |>
  as.data.table()

Study system

plot(arc$timestamp, arc$depth * -1, 
     col = "royalblue", type = "l",
     ylim = c(-max(arc$depth), 0), 
     xlab = "Time stamp", ylab = "Depth (m)")
points(det$timestamp, rep(0, nrow(det)))

Study system

We need to define a map of our study system in R and Julia

Note

We use a coarse raster for convenience only. For a real analysis, a better map would be desirable.

Study system

The map is a raster that sets the spatial domain of our analyses and defines the region(s) within which movements are possible

# Read SpatRaster
map <- dat_gebco()

# dat_gebco() is short-hand for the following:
# spat.tif <- system.file("extdata", "dat_gebco.tif", package = "patter", mustWork = TRUE)
# map      <- terra::rast(spat.tif)

# Export SpatRaster as `env` to Julia
set_map(map)

Study system

Note

In patter v2.0.0, initial particles are sampled in Julia and set_map() includes .as_Raster and .as_GeoArrays arguments if you need to distinguish the two maps. Use .as_Raster = TRUE to define the initial map from which particles are sampled. Use .as_GeoArrays = TRUE to define the map for the movement model.

Study system

Cell entries should be floats
Code inhospitable habitats (e.g., land) as NA_real_
A planar (e.g. UTM) coordinate reference system is required (units: metres)

Note

The planar raster format is used for computational efficiency. In our applications, the map is often a bathymetry SpatRaster. The next version of patter may include routines to facilitate the creation of SpatRaster if you have a shapefile (or similar). For now, see terra::rasterize() to convert a shapefile to a SpatRaster. Raster resolution can be as high as desired, providing the data fits in memory. There is no (minimal?) speed penalty for higher resolution rasters.

Study system

Define the study timeline:
- This sets the period, and resolution, at which we represent animal movement
- We may have observations regularly or irregularly along this timeline

Tip

The appropriate resolution of the timeline is study specific. Considerations include the study species’ movement capacity, the complexity of boundaries (e.g., islands), the resolution at which data were collected and available computational resources. In passive acoustic telemetry studies, a reasonable starting point is to (approximately) align the resolution of the timeline with the transmission interval (≈2 mins).

Study system

Use assemble_timeline() to set a timeline based on the input datasets:

timeline <- assemble_timeline(.datasets = list(det, arc), 
                              .step = "2 mins", 
                              .trim = TRUE)
str(timeline)

 POSIXct[1:24225], format: "2016-03-17 01:50:00" "2016-03-17 01:52:00" "2016-03-17 01:54:00" ...

Study system

Set the timeline using seq.POSIXt():

timeline <- seq(as.POSIXct("2016-03-17 01:50:00", tz = "UTC"),
                as.POSIXct("2016-03-18 01:48:00", tz = "UTC"), 
                by = "2 mins")
str(timeline)

 POSIXct[1:720], format: "2016-03-17 01:50:00" "2016-03-17 01:52:00" "2016-03-17 01:54:00" ...

Tip

Depending on the time step, it is generally wise to split the time series into chunks (e.g., one-month in duration) and run the algorithms for each chunk. Running the algorithms over prolonged periods increases the chances of convergence failure and may exhaust available memory (though patter v2.0.0 can handle this via batching). Here, we use a short, bespoke timeline for example speed only.

Structures

We have defined our study system (map, timeline, observations)
To use this information to reconstruct individual’s movements, our next step is to define the components (structures) of our model:
- State
- Movement model
- Observation model(s)

State

State is an Abstract Type in Patter.jl
This is an object that holds the parameters describing the ‘state’ of an individual
Typically, state means ‘location’ (in two or three dimensions), but additional parameters may be included

# The Julia code for the StateXY struct
struct StateXY <: State
    map_value::Float64
    x::Float64
    y::Float64
end

State

In Patter.jl, the following sub-types are built-in:
- StateXY(map_value, x, y): Used for 2D states (x, y);
- StateXYZ(map_value, x, y, z): Used for 3D states (x, y, z);
- StateCXY(map_value, x, y, heading): Used for correlated 2D states (x, y, heading);
- StateCXYZ(map_value, x, y, z, heading): Used for correlated 3D states (x, y, z, heading);

Tip

It is straightforward new state types; see the help for ?State (or submit a GitHub issue).

State

In patter, you just need to specify the State sub-type as a character string:
- "StateXY" maps onto StateXY in Julia
- "StateXYZ"
- "StateCXY"
- "StateCXYZ"

state <- "StateXY"

Each State sub-type maps onto a corresponding movement model

Movement model

ModelMove is an Abstract Type in Patter.jl
This is an object that holds the dimensions of the movement model, e.g.,
- The map
- A distribution of step lengths
- A distribution of headings

# The Julia code for the ModelMoveXY struct
struct ModelMoveXY{T, U, V, W} <: ModelMove
    map::T
    mobility::U
    dbn_length::V
    dbn_heading::W
end

Movement model

The movement model can be instantiated via an R model_move_*() wrapper
- StateXY maps onto ModelMoveXY, which is wrapped by model_move_xy()
- StateXYZ maps onto ModelMoveXYZ, which is wrapped by model_move_xyz()
- StateCXY maps onto ModelMoveCXY, which is wrapped by model_move_cxy()
- StateCXYZ maps onto ModelMoveCXYZ, which is wrapped by model_move_cxyz()
The fields of the movement model in the model_move_*() wrapper must be specified using Julia code

Movement model

Define a two-dimensional random walk (x, y):

# Define maximum moveable speed between two time steps
mobility   <- 750.0

# Define 2D RW (XY)
model_move <- model_move_xy(.mobility = mobility, 
                            .dbn_length = glue("truncated(Gamma(1.0, 250.0), upper = {mobility})"),
                            .dbn_heading = "Uniform(-pi, pi)")

Tip

In Julia, a number without a decimal point is an integer—which may occasionally cause trouble if a float is required. Specify floats by adding the .0.

Movement model

move_*() wrappers return a string of Julia code:

model_move

ModelMoveXY(env, 750, truncated(Gamma(1.0, 250.0), upper = 750), Uniform(-pi, pi));

Tip

To check the parametrisation of other distributions, such as the Normal distribution in Julia, use JuliaCall::julia_help("Normal"). For new types of movement model; see the help for ?State (or submit a GitHub issue).

Movement model

Visualise the components of our movement model
All built-in movement models have an associated plot method (see ?plot.ModelMove)

plot(model_move)

Movement model

Visualise the components of our movement model manually:

expr <- expression({
  pp <- par(mfrow = c(1, 2))
  curve(dtrunc(x, spec = "gamma", a = 0, b = mobility, shape = 1.0, scale = 250.0), 
        from = 0, to = mobility, n = 1e3L,
        xlab = "Step length (m)", ylab = "Density")
  curve(dunif(x, -pi, pi),
        from = -pi - 0.1, to = pi + 0.1, n = 1e3L,
        xlab = "Heading (rad)", ylab = "Density")
  par(pp)
})

Movement model

Visualise the components of our movement model manually:

eval(expr)

Movement model

Visualise realisations of the movement model:

paths <- sim_path_walk(.map = map, 
                       .timeline = timeline,
                       .state = state,
                       .model_move = model_move, 
                       .n_path = 2L, .one_page = TRUE)

Tip

sim_path_walk() can be used to simulate new trajectories or visualise a prior. The simulation is much faster and more flexible than many alternative routines (e.g., glatos::crw_in_poly()).

Important

In the particle filter, the movement model is a prior. It is a good idea to visualise realisations of the prior. The simulated trajectories should look sensible.

Movement model

Visualise realisations of the movement model:

Movement model

Visualise realisations of the movement model:

head(paths)

   path_id timestep           timestamp map_value        x       y
     <int>    <int>              <POSc>     <num>    <num>   <num>
1:       1        1 2016-03-17 01:50:00  56.13932 698642.1 6267507
2:       1        2 2016-03-17 01:52:00  40.92620 698180.4 6267720
3:       1        3 2016-03-17 01:54:00  45.35426 698061.7 6267769
4:       1        4 2016-03-17 01:56:00  50.36387 698101.0 6267855
5:       1        5 2016-03-17 01:58:00  50.36387 698125.8 6267852
6:       1        6 2016-03-17 02:00:00  40.92620 698116.5 6267734

Tip

Patter.jl automatically truncates the movement model by barriers (e.g., land); that is, we sample movements from a restricted movement model (in the filter) and we evaluate the probability density of truncated moves (in the smoother). With custom movement models, you need to specify Patter.simulate_step() method (for the filter) and a Patter.logpdf_step() (for the smoother). These methods simulate and evaluate the log-probability of an unrestricted step. Patter.jl handles the truncation for you (by Patter.simulate_move() and Patter.logpdf_move()).

Movement model

More sophisticated movement models are possible, including:
- Random walks in two or three dimensions
- Correlated random walks
- Behavioural-switching models
See the help for ?State to get started

Movement model

This is an example of a more sophisticated movement model for flapper skate:

Example behavioural-switching model

Movement model

There is a trade-off between simplicity and complexity
It is convenient to assume that static parameters are known
Where parameters are unknown, you can:
1. See Lavender et al. (2025a) for guidance on parameter sensitivity
2. Run a sensitivity analysis:
  - Use simulations
  - Compare results from real-world analyses based on a small number of sensible models
3. Formally propagate uncertainty with multiple algorithm runs
4. Data-driven parameter estimation (if data are sufficient)

Observation model(s)

We have defined the movement model
This encodes the process of animal movement
We now need to define the observation models
These encode the probability of the observations, given the latent (unobserved) location of the animal

Observation model(s)

As a reminder, we have collected detections and depth time series:

head(det, 3)

             timestamp receiver_id
                <POSc>       <int>
1: 2016-03-17 01:50:00          26
2: 2016-03-17 01:52:00          26
3: 2016-03-17 01:54:00          26

head(arc, 3)

             timestamp  depth
                <POSc>  <num>
1: 2016-03-16 00:00:00 159.15
2: 2016-03-16 00:02:00 161.01
3: 2016-03-16 00:04:00 162.86

Observation model(s)

Each dataset should comprise a time series of observations for one individual
Datasets should be pre-processed as required: e.g.,
- Check receiver locations are valid on the map
- Identify false detections
patter accepts any kind of observational dataset

Tip

For passive acoustic telemetry and/or archival data, the pat_setup_data() function is available to validate the format of these data types for use in patter.

Observation model(s)

Each dataset should be provided as a data.table with the following columns:
- sensor_id, timestamp, obs;
- Additional columns that define observation-model parameters;

Tip

See ?assemble for full details. Note that the resolution of the time stamps must match the resolution of the timeline, although assemble_custom() can handle this for you (but irregular observations are fine).

Observation model(s)

In passive acoustic telemetry datasets, we record detections but the observations comprise detections and non-detections
Use the assemble_acoustics() helper to assemble the acoustic (detection, non-detection) time series:

acc <- assemble_acoustics(.timeline   = timeline,
                          .detections = det,
                          .moorings   = dat_moorings, 
                          .services   = NULL)

Observation model(s)

head(acc)

             timestamp sensor_id   obs receiver_x receiver_y receiver_alpha
                <POSc>     <int> <int>      <num>      <num>          <num>
1: 2016-03-17 01:50:00         3     0   706442.1    6254007              4
2: 2016-03-17 01:50:00         4     0   709742.1    6267707              4
3: 2016-03-17 01:50:00         7     0   708742.1    6269107              4
4: 2016-03-17 01:50:00         9     0   706042.1    6254307              4
5: 2016-03-17 01:50:00        11     0   707542.1    6267707              4
6: 2016-03-17 01:50:00        12     0   710042.1    6267307              4
   receiver_beta receiver_gamma
           <num>          <num>
1:         -0.01            750
2:         -0.01            750
3:         -0.01            750
4:         -0.01            750
5:         -0.01            750
6:         -0.01            750

Note

assemble_acoustics() accounts for variable receiver deployment periods and servicing events. Duplicate observations (that is, detections at the same receiver in the same time step) are dropped. If available in .moorings, additional columns (receiver_alpha, receiver_beta and receiver_gamma) are included as required for the default acoustic observation model (ModelObsAcousticLogisTrunc). But you have the flexibility to define receiver- and time-varying parameters by modifying the data.table as required.

Observation model(s)

If you have passive acoustic telemetry observations, you should also use acoustic containers
The syntax is:

# ModelObsContainer
containers <- assemble_acoustics_containers(.timeline, .acoustics, .mobility, ...)

Note

Acoustic containers are a special kind of ‘observation’ that enhance the efficiency of the particle filter (reducing the number of particles required to achieve convergence). For simplicity, we don’t implement them here.

Observation model(s)

We have to assemble other datasets with the required columns (sensor_id, timestamp, obs, ...) manually

head(arc, 2)

             timestamp  depth
                <POSc>  <num>
1: 2016-03-16 00:00:00 159.15
2: 2016-03-16 00:02:00 161.01

arc <- 
  arc |> 
  # Define 'dummy' sensor_id
  mutate(sensor_id = 1L) |> 
  # Observations must be in the 'obs' column
  rename(obs = "depth") |> 
  # Select columns
  select("sensor_id", "timestamp", "obs") |>
  # Add observation model parameters 
  mutate(depth_sigma = 50, depth_deep_eps = 20) |>
  as.data.table()

Observation model(s)

It is a good idea to pass other datasets through the assemble_custom() function for processing/checks though

arc <- assemble_custom(.timeline = timeline, .dataset = arc)
head(arc)

   sensor_id           timestamp   obs depth_sigma depth_deep_eps
       <int>              <POSc> <num>       <num>          <num>
1:         1 2016-03-17 01:50:00 73.78          50             20
2:         1 2016-03-17 01:52:00 73.32          50             20
3:         1 2016-03-17 01:54:00 73.32          50             20
4:         1 2016-03-17 01:56:00 73.32          50             20
5:         1 2016-03-17 01:58:00 73.55          50             20
6:         1 2016-03-17 02:00:00 68.70          50             20

Observation model(s)

For each dataset, we need to define the corresponding ModelObs sub-type
A ModelObs sub-type is a structure that holds the parameters of an observation model
ModelObs sub-types are instantiated in Julia using the data.table(s) assembled above
We use ModelObs instances to simulate observations or to evaluate the log-probability of observations, given the latent location

Observation model(s)

The following ModelObs structures are built-in to Patter.jl:
- ModelObsAcousticLogisTrunc
- ModelObsDepthUniformSeabed
- ModelObsDepthNormalTruncSeabed
- ModelObsContainer

Tip

It is straightforward to define new ModelObs sub-types for new data types; see ?ModelObs (or submit a GitHub issue).

Observation model(s)

ModelObsAcousticLogisTrunc is designed for acoustic observations
The structure holds the parameters of truncated logistic detection probability model
This contains the fields we need to simulate an observation, or evaluate the log-probability of an observation, given a transmission from a particular location

Observation model(s)

# Julia code for the ModelObsAcousticLogisTrunc sub-type
struct ModelObsAcousticLogisTrunc <: ModelObs
    sensor_id::Int64
    receiver_x::Float64
    receiver_y::Float64
    receiver_alpha::Float64
    receiver_beta::Float64    
    receiver_gamma::Float64
end

Observation model(s)

In a truncated logistic detection probability model:
- receiver_alpha and receiver_beta affect the rate of decline in detection probability with distance from a receiver
- receiver_gamma is the maximum detection range

Observation model(s)

Example truncated logistic detection probability models

Observation model(s)

You can explore different parameters values using this code:

expr <- expression({
  # Define parameters 
  receiver_alpha <- 4
  receiver_beta  <- -0.01
  receiver_gamma <- 1000
  # Visualise model 
  curve(ifelse(x <= receiver_gamma, 
              1 / (1 + exp(-(receiver_alpha + receiver_beta * x))), 
              0), 
       from = 0, to = receiver_gamma, 
       xlab = "Distance (m)", ylab = "Detection probability")
})

Observation model(s)

eval(expr)

Observation model(s)

ModelObsDepthNormalTruncSeabed is a ModelObs structure for a depth observation and a truncated normal model
This is a suitable depth-observation model for benthic/demersal species
The individual must be located in an envelope around the bathymetric depth, defined by a normal distribution centred at this location, truncated at depth_deep_eps m below the seabed and 0 m, with standard deviation depth_sigma

Observation model(s)

Example truncated normal depth observation models

Observation model(s)

ModelObsDepthUniformSeabed is ModelObs structure for a depth observation and a uniform depth model
This model assumes that an individual must be located in an envelope around the bathymetric depth, defined by two error terms (depth_shallow_eps and depth_deep_eps)
This model can be tailored for benthic or pelagic species

Observation model(s)

ModelObsContainer is a ModelObs structure for capture/recapture/acoustic containers
This is a special ‘observation’ type that represents the distance the individual must be from some future location (e.g., a recapture location)
ModelObsContainers are effectively a redundant use of the data that mitigates particle degeneracy

Observation model(s)

For the particle filter, multiple sets of observations should be provided as a named list
List names define the model sub-types

yobs <- list(ModelObsAcousticLogisTrunc = acc, 
             ModelObsDepthNormalTruncSeabed = arc)

Note

In Julia, the list of observations is translated into a hash table of ModelObs sub-types.

Observation model(s)

As for the movement model, the parameters of the observation model(s) are specified by the user
See the earlier advice for situations where parameters are uncertain

Particle filter

We have now defined:
1. Study area map
2. Study timeline
3. State
4. Movement model
5. Observations/observation model parameters/observation model sub-types
We are now in a position to run the particle filter!

Particle filter

At $t = 1$:
1. Sample initial particles (locations → States), via Patter.simulate_states_init()
2. For each simulated particle, compute the (log) weight, via Patter.logpdf_obs()
3. Record a subset of .n_record resampled (equally weighted) particles in memory

Particle filter

For $t \in 2, ..., T$
1. Simulate animal movement, generating new particles (locations → States), via Patter.simulate_step()
2. For each simulated particle, compute the (log) weight, via Patter.logpdf_obs()
3. Record a subset of .n_record resampled (equally weighted) particles in memory
4. Periodically, resample particles with replacement (using the low-variance resampling algorithm)

Tip

For custom State, ModelMove and/or ModelObs structures, associated methods for corresponding Julia routines are required (see ?State, ?ModelMove and ?ModelObs for guidance).

Particle filter

The particle filter is implemented via pf_filter() with the following arguments:

pf_filter(
  .timeline,
  .state = "StateXY",
  .xinit = NULL,
  .model_move = model_move_xy(),
  .yobs,
  .n_move = 100000L,
  .n_particle = 1000L,
  .n_resample = as.numeric(.n_particle),
  .t_resample = NULL,
  .n_record = 1000L,
  .n_iter = 1L,
  .direction = c("forward", "backward"),
  .batch = NULL,
  .collect = TRUE,
  .progress = julia_progress(),
  .verbose = getOption("patter.verbose")
)

Particle filter

Argument	Description
`.timeline`	The timeline
`.state`	The state
`.xinit`	Initial location(s)
`.model_move`	Movement model
`.yobs`	The named list of observations
`.n_move`	Set to default
`.n_particle`	The number of particles
`.n_resample`	Set to default
`.t_resample`	Set to default
`.n_record`	Set to default
`.n_iter`	Number of particle filter iterations
`.direction`	The direction of the particle filter
`.batch`	Implement batching to handle memory requirements
`.collect`	Collect outputs from `Julia` in `R`
`.progress`	Progress bar options
`.verbose`	Show user output

Particle filter

To set initial states:
1. Specify known states via .xinit (or via set_map()’s .as_Raster argument)
2. Let patter sample states automatically from part of the map that is compatible with the initial observations

Tip

If you know the initial location(s) of the animal, provide this information via .xinit (or via set_map()’s .as_Raster argument). If .xinit is not provided, initial locations are sampled uniformly from the part of the map that is compatible with the initial observations. In patter v2.0.0, this happens in Julia. You can set a unique initial map via set_map(..., .as_Raster = TRUE, .as_GeoArray = FALSE). For custom ModelObs subtypes, Patter.map_init() methods are required for automated sampling of the part of the map compatible with observations. In patter v2.0.0, methods written in Julia are required.

Particle filter

Let’s set up the filter for our example time series:

# List filter arguments
args <- list(.timeline    = timeline,
             .state       = state,
             .yobs        = yobs,
             .model_move  = model_move,
             .n_particle  = 1e5L)

Tip

How many particles do we need? Recall that the particles approximate the distribution of an individual’s State at every time step ($f(\mathbf{s}_t \mid \mathbf{y}_{1:t})$ in the filter). The number of particles required to approximate a distribution depends on its dimensionality. For a two-dimensional distribution, ≈1000 particles may be sufficient. In practice, in the filter, a larger number of particles may be required to ensure that enough particles remain ‘alive’ at each time step to approximate the distribution. With acoustic observations alone, relatively few particles are required. (For passive acoustic telemetry data, the ModelObsContainer structure facilitates convergence with fewer particles than used here (≤ 25,000)). With additional observations, more particles may be required depending on the complexity of the landscape that particles have to explore. This is associated with a linear speed penalty.

Particle filter

Let’s run the filter!

# Forward run
fwd <- do.call(pf_filter, args, quote = TRUE)

Warning

We are approximating the (partial) marginal distribution of the individual’s location, given the data up to and including that time, at every time step (every two minutes here!). This is awesome and the Julia code is very fast, but it is a lot to ask and for large numbers of particles, observations and time steps, the filter may take minutes or hours to run (especially on older computers).

Particle filter

The filter returns a pf_particles-class object
This is a named list:

  class(fwd)

[1] "list"         "pf_particles"

  summary(fwd)

            Length Class      Mode
states      6      data.table list
diagnostics 4      data.table list
callstats   7      data.table list

Particle filter

states is a data.table of particles:

head(arrange(fwd$states, timestep))

   path_id timestep           timestamp map_value        x       y
     <int>    <int>              <POSc>     <num>    <num>   <num>
1:       1        1 2016-03-17 01:50:00  81.56814 709342.1 6253607
2:       2        1 2016-03-17 01:50:00  58.17422 709242.1 6253107
3:       3        1 2016-03-17 01:50:00  87.54828 709142.1 6253407
4:       4        1 2016-03-17 01:50:00 113.60658 708742.1 6253207
5:       5        1 2016-03-17 01:50:00  83.91008 708942.1 6253107
6:       6        1 2016-03-17 01:50:00  69.39536 709442.1 6253607

Warning

The path_id is simply an index; it is named path_id for consistency with sim_path_walk(). These are not trajectories.

Particle filter

Map particle coordinates:

plot_xyt(.map = map, .coord = fwd$states, .steps = 1L, 
         .add_points = list(pch = ".", col = "red"))

Tip

In pf_filter(), at each time step, .n_record particles are sampled with replacement, in line with the weights, and saved in memory. All recorded particles thus have equal weight.

Tip

See the help for ?plot_xyt() to create animations.

Particle filter

We can visualise the probability distribution for the individual’s initial location as the proportion of particles in each grid cell:

map_pou(.map = map, .coord = fwd$states[timestep == 1L, ])

Note

Recall that particles are equally weighted at this point thanks to resampling.

Tip

map_*() functions such as map_pou() and map_dens() accept any data.table of (x, y) coordinates. Time series and weights (marks) are handled automatically.

Particle filter

diagnostics is a data.table of filter diagnostics:

fwd$diagnostics

     timestep           timestamp      ess     maxlp
        <int>              <POSc>    <num>     <num>
  1:        1 2016-03-17 01:50:00 42322.37 -4.640798
  2:        2 2016-03-17 01:52:00 55025.08 -4.620056
  3:        3 2016-03-17 01:54:00 60575.36 -4.619843
  4:        4 2016-03-17 01:56:00 49133.86 -4.242834
  5:        5 2016-03-17 01:58:00 46250.73 -4.621105
 ---                                                
716:      716 2016-03-18 01:40:00 80784.34 -4.408220
717:      717 2016-03-18 01:42:00 85413.20 -4.408251
718:      718 2016-03-18 01:44:00 88642.58 -4.408215
719:      719 2016-03-18 01:46:00 88424.21 -4.408220
720:      720 2016-03-18 01:48:00 88707.48 -4.408224

Particle filter

maxlp is the maximum log-posterior at each time step (assuming resampling @ every time step):

plot(fwd$diagnostics$timestamp, fwd$diagnostics$maxlp, 
     xlab = "Time stamp", ylab = "maxlp",
     type = "l")

Tip

exp(maxlp) is the highest likelihood score at each time step (0–1). Watch out for time steps with very low maxlp. For passive acoustic telemetry analyses, ModelObsAcousticLogisTrunc helps to prevent this via the receiver_gamma threshold.

Particle filter

ess is the effective sample size (ESS):

plot(fwd$diagnostics$timestamp, fwd$diagnostics$ess, 
     xlab = "Time stamp", ylab = "ESS",
     type = "l")
points(det$timestamp, rep(0, nrow(det)),
       pch = 21, col = "red", bg = "red", cex = 0.5)
abline(h = 500, col = "royalblue", lty = 3)

Tip

For a 2D distribution (which we have here), ideally the ESS should be >= 500 particles to achieve a reasonable approximation. Small ESS may suggest we need to increase the number of particles (or revise our assumptions). However, small ESS is not necessarily problematic, as there may be time steps when there were only a small number of possible locations in which the individual could have been located.

Particle filter

callstats defines call statistics
This includes a column for log-likelihood of the filter run (loglik) and whether or not the filter converged (convergence):

fwd$callstats

             timestamp         routine n_particle n_iter    loglik convergence
                <POSc>          <char>      <int>  <int>     <num>      <lgcl>
1: 2025-06-01 13:50:17 filter: forward     100000      1 -3409.552        TRUE
       time
      <num>
1: 7.858565

Warning

Convergence failures produce a warning. These may occur for a variety of reasons including (1) code bugs, (2) prior–data conflicts, (3) incorrectly specified observation models and (4) extreme particle degeneracy. With acoustic observations alone, likely causes are 1–3. We have observed extreme particle degeneracy (4) in complex, labyrinthine landscapes, where particles have to explore a large number of possible routes of which only a tiny fraction are ultimately compatible with data.

Particle filter

The particle filter approximates the partial marginal distribution of the individual’s location given the data up to and including each time step, that is, $f(\mathbf{s}_t \mid \mathbf{y}_{1:t})$
If we run the particle filter forwards and backwards, we can use the two-filter smoother to account for all data, that is, $f(\mathbf{s}_t \mid \mathbf{y}_{1:T})$
In most realistic scenarios, this substantially improves maps of space use

Particle filter

Run the filter backwards:

# Backward run
args$.direction <- "backward"
bwd <- do.call(pf_filter, args, quote = TRUE)

Tip

Take care when specifying initial locations for the backward run.

Two-filter smoother

Smoothing requires us to compute the probability density of movements between locations
We account for automatic truncation of the movement model via a Monte Carlo simulation
To speed up the simulation, we can set a ‘validity map’ for certain kinds of movement model
This is the region within which movements are always valid

Two-filter smoother

pp <- par(mfrow = c(1, 2))
terra::plot(map)
set_vmap(.map = map, .mobility = mobility, .plot = TRUE)

par(pp)

Two-filter smoother

Run the two-filter smoother using outputs from pf_filter() in the Julia environment:

# Run smoother 
smo <- pf_smoother_two_filter(.n_particle = 100L, .n_sim = 100L)

# Reset vmap
set_vmap()

Two-filter smoother

Argument	Description
`.n_particle`	The number of particles from the filter used in the algorithm
`.n_sim`	The number of Monte Carlo simulations used to compute the normalisation constant (when required)
`.cache`	Cache normalisation constants for speed
`.batch`	Use batching to handle large memory requirements
`.collect`	Collect outputs in Julia in R
`.progress`	Progress bar options
`.verbose`	Show user output

Two-filter smoother

Note

Smoothing is expensive ($\mathcal{O}(TN^2)$). In this example, we use .n_particle = 100L for speed only. For real analyses, .n_particle = 1000L is generally acceptable. .n_sim = 100L is reasonable.

Note

Smoothing uses the probability density of movements between pairs of locations. The density of an unrestricted move is evaluated by Patter.logpdf_step(). For custom movement models, a Patter.logpdf_step() method is required.

Two-filter smoother

pf_smoother_two_filter() returns the familiar pf_particles-class object:

class(smo)

[1] "list"         "pf_particles"

summary(smo)

            Length Class      Mode
states      6      data.table list
diagnostics 4      data.table list
callstats   7      data.table list

Two-filter smoother

We can analyse the outputs in the same way as those from the filter:

head(arrange(smo$states, timestep))

   path_id timestep           timestamp map_value        x       y
     <int>    <int>              <POSc>     <num>    <num>   <num>
1:       1        1 2016-03-17 01:50:00 121.18092 708843.5 6253359
2:       2        1 2016-03-17 01:50:00  93.25472 708905.1 6253228
3:       3        1 2016-03-17 01:50:00  74.88351 709052.7 6253112
4:       4        1 2016-03-17 01:50:00  99.90396 709006.6 6253401
5:       5        1 2016-03-17 01:50:00  83.21526 709091.9 6253249
6:       6        1 2016-03-17 01:50:00  80.03239 709163.3 6253289

Two-filter smoother

diagnostics includes ESS:

head(smo$diagnostics)

   timestep           timestamp       ess maxlp
      <int>              <POSc>     <num> <num>
1:        1 2016-03-17 01:50:00 100.00000   NaN
2:        2 2016-03-17 01:52:00  84.46088   NaN
3:        3 2016-03-17 01:54:00  82.85654   NaN
4:        4 2016-03-17 01:56:00  54.97102   NaN
5:        5 2016-03-17 01:58:00  70.09009   NaN
6:        6 2016-03-17 02:00:00  73.17920   NaN

Two-filter smoother

callstats:

smo$callstats

             timestamp              routine n_particle n_iter loglik
                <POSc>               <char>      <int>  <int>  <num>
1: 2025-06-01 13:50:33 smoother: two-filter        100     NA    NaN
   convergence     time
        <lgcl>    <num>
1:       FALSE 0.795193

Mapping

We can compute probability-of-use (POU) using particle samples from the filter (for comparison) or the smoother (recommended):

map_pou(.map = map, .coord = smo$states)$ud

Mapping

In practice, POU can be sensitive to grid resolution and kernel smoothing may be desirable

map_dens(
  .map,
  .owin = as.owin.SpatRaster(.map),
  .coord = NULL,
  .discretise = FALSE,
  .shortcut = list(),
  .sigma = bw.h,
  ...,
  .fterra = FALSE,
  .plot = TRUE,
  .tryCatch = TRUE,
  .verbose = getOption("patter.verbose")
)

Mapping

map_dens() implements kernel smoothing via spatstat.explore::density.ppp()
- Use .sigma to choose the method to estimate the bandwidth
- E.g., bw.diggle implements cross validation (expensive)
- Edge-correction is automatically implemented

Note

The default sigma in map_dens() differs from density.ppp(). map_dens() uses bw.h(), which sets the bandwidth based on the combined variance of summarised coordinates using the ad-hoc method (Worton, 1989); i.e., bw.h = function(X) sqrt(0.5 * (var(X$x) + var(X$y))) * (X$n^(-1/6)).

Tip

To improve speed: specify .owin, using a simplified sf window, set .discretise = TRUE and use .fterra = TRUE.

Mapping

Implement kernel smoothing with smoothed particles:

# Estimate UD
ud <- map_dens(.map = map,
               .coord = smo$states,
               .sigma = bw.diggle)$ud

# Add home range
map_hr_home(ud, .add = TRUE)

Important

This map is a biologically and probabilistically sound utilisation distribution. Cell colours denote the probability density of observing an individual in a given grid cell at a randomly chosen time.

Warning

Don’t confuse particle smoothing with post-hoc kernel smoothing!

Mapping

Residency

We can easily compute residency in a particular area (at any arbitrary time-scale of our choosing) as the proportion of particles inside/outside an area:

# Define an example polygon (e.g., MPA)
spatpoly <- 
  cbind(708433.9, 6251765) |> 
  terra::vect(crs = terra::crs(map)) |> 
  terra::buffer(width = 2000)

# Visualise study area and 'MPA'
terra::plot(map)
points(smo$states$x, smo$states$y, pch = ".", col = "red")
terra::lines(spatpoly, col = "black", lwd = 2)

Residency

# Lookup whether or not particles are in the polygon
spatcoord <- terra::vect(cbind(smo$states$x, smo$states$y))
smo$states[, in_poly := 
               terra::relate(spatcoord, spatpoly, relation = "intersects")]
               
# Calculate the expected time spent in the 'MPA' (%)
table(smo$states$in_poly) / nrow(smo$states) * 100


  FALSE    TRUE 
71.8875 28.1125

Important

This estimate of residency is probabilistically sound.

Recap

Boilerplate setup (1/7)

# Load & attach package 
library(patter)

# Load datasets
head(det, 2)
head(arc, 2)

# Connect to Julia 
julia_connect()
set_seed()

# Set map 
map <- dat_gebco()
set_map(map)

# Study timeline 
timeline <- assemble_timeline(.datasets = list(det, arc))

Set state (2/7)

state <- "StateXY"

Set movement model (3/7)

mobility   <- 750.0
model_move <- model_move_xy(.mobility = mobility, 
                            .dbn_length = glue("truncated(Gamma(1.0, 250.0), upper = {mobility})"),
                            .dbn_heading = "Uniform(-pi, pi)")

Set observation model(s) (4/7)

# Acoustic observations 
acc <- assemble_acoustics(.timeline   = timeline,
                          .detections = det,
                          .moorings   = dat_moorings, 
                          .services   = NULL)

# Archival observations 
arc <- assemble_custom(.timeline = timeline, 
                        .archival = 
                           arc |> 
                             rename(obs = "depth") |> 
                             select("sensor_id", "timestamp", "obs") |>
                             mutate(depth_sigma = 50, depth_deep_eps = 20) |> 
                           as.data.table())

# All observations
yobs <- list(ModelObsAcousticLogisTrunc = acc, 
             ModelObsDepthNormalTruncSeabed = arc)

Run particle filter (5/7)

# List filter arguments
args <- list(.map         = map,
             .timeline    = timeline,
             .state       = state,
             .yobs        = yobs,
             .model_move  = model_move,
             .n_particle  = 1e5L)

# Forward run
fwd <- do.call(pf_filter, args, quote = TRUE)

# Backward run
args$.direction <- "backward"
bwd <- do.call(pf_filter, args, quote = TRUE)

Run particle smoother (6/7)

smo <- pf_smoother_two_filter()

Analyse outputs (7/7)

map_dens(.map = map, .coord = smo$states, .sigma =  bw.diggle)

Simulation-based workflow

patter also supports simulation-based workflows:
1. sim_array() simulates arrays
2. sim_path_walk() simulates movement paths
3. sim_observations() simulates observations
You can model simulated observations

Important

Simulation studies inform real-world analyses.

Simulate arrays

Use sim_array() to simulate randomly or regularly arranged receivers
Locations are sampled from .map
Default observation model parameters (receiver_alpha, receiver_beta, receiver_gamma) can be assigned

moorings <- sim_array(.map = map, 
                      .timeline = timeline,
                      .arrangement = "random",
                      .n_receiver = 100L, 
                      .receiver_alpha = 5, 
                      .receiver_beta = -0.01, 
                      .receiver_gamma = 1200, 
                      .n_array = 1L)

Simulate arrays

head(moorings)

   array_id receiver_id      receiver_start        receiver_end receiver_x
      <int>       <int>              <POSc>              <POSc>      <num>
1:        1           1 2016-03-17 01:50:00 2016-03-18 01:48:00   708842.1
2:        1           2 2016-03-17 01:50:00 2016-03-18 01:48:00   706142.1
3:        1           3 2016-03-17 01:50:00 2016-03-18 01:48:00   708942.1
4:        1           4 2016-03-17 01:50:00 2016-03-18 01:48:00   705542.1
5:        1           5 2016-03-17 01:50:00 2016-03-18 01:48:00   704542.1
6:        1           6 2016-03-17 01:50:00 2016-03-18 01:48:00   709542.1
   receiver_y receiver_alpha receiver_beta receiver_gamma
        <num>          <num>         <num>          <num>
1:    6254807              5         -0.01           1200
2:    6253307              5         -0.01           1200
3:    6253107              5         -0.01           1200
4:    6261407              5         -0.01           1200
5:    6252907              5         -0.01           1200
6:    6268407              5         -0.01           1200

Simulate paths

We have already seen how to simulate paths

# Set state and movement model 
state      <- "StateXY"
mobility   <- 500.0
model_move <- 
  model_move_xy(.mobility = mobility, 
                .dbn_length = glue("truncated(Normal(50, 250.0), lower = 0, upper = {mobility})"), 
                .dbn_heading = "VonMises(0, 0.25)")

# Simulate path using movement model 
pp <- par(mfrow = c(1, 2))
path_coord <- sim_path_walk(.map = map, 
                            .timeline = timeline, 
                            .state = state, 
                            .n_path = 1L, 
                            .model_move = model_move)

# Fit UD using cross validation
path_ud <- map_dens(.map = map, 
                    .coord = path_coord, 
                    .discretise = TRUE,
                    .fterra = TRUE, 
                    .sigma =  bw.diggle)$ud
par(pp)

Simulate paths

head(path_coord)

   path_id timestep           timestamp map_value        x       y
     <int>    <int>              <POSc>     <num>    <num>   <num>
1:       1        1 2016-03-17 01:50:00  19.86056 704742.1 6262707
2:       1        2 2016-03-17 01:52:00  23.43867 704559.9 6263085
3:       1        3 2016-03-17 01:54:00  23.43867 704540.4 6263136
4:       1        4 2016-03-17 01:56:00  29.62455 704605.7 6263275
5:       1        5 2016-03-17 01:58:00  38.94438 704805.0 6263067
6:       1        6 2016-03-17 02:00:00  26.42981 704543.1 6263332

Simulate observations

To simulate observations, we need to choose our observation models and their parameters
For each observation model, we need to provide parameters as a data.table that defines, for each sensor_id (e.g., receiver), the model parameters associated with that sensor

Note

Time-varying observation model parameters are not currently supported in sim_observations() (but are fully supported in pf_filter()).

Simulate observations

Let’s set things up to simulate acoustic observations:

# ModelObsAcousticLogisTrunc 
# > We have multiple acoustic sensors
# > This data.table includes the relevant parameters
pars_ModelObsAcousticLogisTrunc <-
    moorings |>
    select(sensor_id = "receiver_id",
           "receiver_x", "receiver_y",
           "receiver_alpha", "receiver_beta", "receiver_gamma") |>
    as.data.table()

Simulate observations

And depth observations:

# ModelObsDepthNormalTruncSeabed
# > We have one depth sensor 
# > This data.table includes the relevant parameters 
pars_ModelObsDepthNormalTruncSeabed <- data.frame(sensor_id = 1L,
                                                  depth_sigma = 20,
                                                  depth_deep_eps = 20)

Tip

Use ?ModelObs or JuliaCall::julia_help("ModelObs") for information on built-in ModelObs structures and their parameters.

Tip

Of course, we can define and use custom ModelObs sub-types. Just define the sub-type and an associated Patter.simulate_obs() method. See the help for ?ModelObs for examples.

Simulate observations

Simulate observations along .timeline via sim_observations():

obs <- sim_observations(
  .timeline  = timeline, 
  .model_obs = list(ModelObsAcousticLogisTrunc = pars_ModelObsAcousticLogisTrunc, 
                    ModelObsDepthNormalTruncSeabed = pars_ModelObsDepthNormalTruncSeabed))

Note

For each sensor, an observation is simulated for each time step. Irregular observations are not not currently supported in sim_observations() (but are fully supported in pf_filter()).

Simulate observations

sim_observations() returns a list, with one element for every .model_obs
Each element is a list, with one sub-element for each simulated path
Each sub-element is a data.table that contains the observations

Simulate observations

str(obs)

List of 2
 $ ModelObsAcousticLogisTrunc    :List of 1
  ..$ :Classes 'data.table' and 'data.frame':   72000 obs. of  8 variables:
  .. ..$ timestamp     : POSIXct[1:72000], format: "2016-03-17 01:50:00" "2016-03-17 01:50:00" ...
  .. ..$ obs           : int [1:72000] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..$ sensor_id     : int [1:72000] 1 2 3 4 5 6 7 8 9 10 ...
  .. ..$ receiver_x    : num [1:72000] 708842 706142 708942 705542 704542 ...
  .. ..$ receiver_y    : num [1:72000] 6254807 6253307 6253107 6261407 6252907 ...
  .. ..$ receiver_alpha: num [1:72000] 5 5 5 5 5 5 5 5 5 5 ...
  .. ..$ receiver_beta : num [1:72000] -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 ...
  .. ..$ receiver_gamma: num [1:72000] 1200 1200 1200 1200 1200 1200 1200 1200 1200 1200 ...
  .. ..- attr(*, ".internal.selfref")=<externalptr> 
 $ ModelObsDepthNormalTruncSeabed:List of 1
  ..$ :Classes 'data.table' and 'data.frame':   720 obs. of  5 variables:
  .. ..$ timestamp     : POSIXct[1:720], format: "2016-03-17 01:50:00" "2016-03-17 01:52:00" ...
  .. ..$ obs           : num [1:720] 12.6 42.1 27 12.8 31 ...
  .. ..$ sensor_id     : int [1:720] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ depth_sigma   : num [1:720] 20 20 20 20 20 20 20 20 20 20 ...
  .. ..$ depth_deep_eps: num [1:720] 20 20 20 20 20 20 20 20 20 20 ...
  .. ..- attr(*, ".internal.selfref")=<externalptr>

Simulate observations

acc <- obs$ModelObsAcousticLogisTrunc[[1]]
head(acc)

             timestamp   obs sensor_id receiver_x receiver_y receiver_alpha
                <POSc> <int>     <int>      <num>      <num>          <num>
1: 2016-03-17 01:50:00     0         1   708842.1    6254807              5
2: 2016-03-17 01:50:00     0         2   706142.1    6253307              5
3: 2016-03-17 01:50:00     0         3   708942.1    6253107              5
4: 2016-03-17 01:50:00     0         4   705542.1    6261407              5
5: 2016-03-17 01:50:00     0         5   704542.1    6252907              5
6: 2016-03-17 01:50:00     0         6   709542.1    6268407              5
   receiver_beta receiver_gamma
           <num>          <num>
1:         -0.01           1200
2:         -0.01           1200
3:         -0.01           1200
4:         -0.01           1200
5:         -0.01           1200
6:         -0.01           1200

Simulate observations

arc <- obs$ModelObsDepthNormalTruncSeabed[[1]]
head(arc)

             timestamp      obs sensor_id depth_sigma depth_deep_eps
                <POSc>    <num>     <int>       <num>          <num>
1: 2016-03-17 01:50:00 12.64076         1          20             20
2: 2016-03-17 01:52:00 42.14758         1          20             20
3: 2016-03-17 01:54:00 27.04716         1          20             20
4: 2016-03-17 01:56:00 12.76650         1          20             20
5: 2016-03-17 01:58:00 31.00878         1          20             20
6: 2016-03-17 02:00:00 35.34954         1          20             20

Simulate observations

# Plot depths
  plot(arc$timestamp, arc$obs * -1, 
       col = "royalblue", type = "l",
       ylim = c(-max(arc$obs), 0), 
       xlab = "Time stamp", ylab = "Depth (m)")
# Add detections 
det <- acc[obs == 1L, ]
points(det$timestamp, rep(0, nrow(det)), col = "red")

Simulate observations

Run algorithms

Now we can run the algorithms and compare the simulated pattern of space use to the reconstructed pattern:

yobs <- list(ModelObsAcousticLogisTrunc = acc,
             ModelObsDepthNormalTruncSeabed = arc)

# List filter arguments
args <- list(.timeline    = timeline,
             .state       = state,
             .yobs        = yobs,
             .model_move  = model_move,
             .n_particle  = 5e4L, 
             .verbose     = FALSE)

# Forward run
fwd <- do.call(pf_filter, args, quote = TRUE)

# Backward run
args$.direction <- "backward"
bwd <- do.call(pf_filter, args, quote = TRUE)

# Run smoother 
smo <- pf_smoother_two_filter(.n_particle = 100)

Mapping

pp <- par(mfrow = c(1, 3))

# Plot the UD for the simulated path 
terra::plot(path_ud)
map_hr_home(path_ud, .add = TRUE)

# Plot the estimated UD 
patter_ud <- map_dens(.map = map, 
                      .coord = smo$states, 
                      .discretise = TRUE, 
                      .fterra = TRUE,
                      .sigma =  bw.diggle)$ud
map_hr_home(patter_ud, .add = TRUE)

# Compare to a UD based on the mean-position (COA) algorithm
coa_ud <- map_dens(.map = map, 
                   .coord = coa(.map          = map, 
                                .detections   = det, 
                                .moorings     = moorings, 
                                .delta_t      = "2 hours", 
                                .plot_weights = FALSE), 
                   .discretise = FALSE, 
                   .fterra = TRUE, 
                   .sigma =  bw.diggle)$ud

par(pp)

Mapping

Ancillary functions

See ?patter for a full list of patter functions
- ?datasets-mefs for example datasets
- ?datasets-algorithms for example outputs
- ?patter-progress for progress options
- ?glossary for a glossary
Useful ancillary functions include:
- pat_setup_data() validates passive acoustic telemetry datasets
- coa() runs the COA algorithm
- map_hr_*() functions estimate home ranges
- skill_*() function evaluate model skill
- cl_lapply() handles efficient R-side parallelisation
- file_*() wrappers interact with the file system

Advanced

Custom `State`, `ModelMove` `ModelObs`

To incorporate custom structures, we need to:
1. Register the structures in Julia
2. Write associated Julia methods
3. Write any R wrappers for the routines

# Julia code
# * Load additional packages for these examples
using Dates

Tip

It is a good idea to include Julia code in ./Julia/src in a patter-structs.jl file (or similar). Source Julia code with JuliaCall::julia_source().

Custom `State` and `ModelMove`

For custom State and ModelMove structures, we need to:
1. Register the structures in Julia
2. (optional) Write a Patter.states_init() methods to translate initial coordinates to States
3. Write simulation and density methods:
  - Patter.simulate_step(state, model, ...) (for simulation)
  - Patter.logpdf_step(state_from, state_to, model_move, ...) (for smoothing)
4. (optional) Write a model_move_*() wrapper for convenience

Note

Julia uses multiple dispatch; that is, the generic simulate_*() and logpdf_*() functions ‘call’ different ‘methods’ depending on the state and model that is provided. This is similar to S4 dispatch in R.

Custom `ModelObs`

For custom ModelObs structures, we need to:
1. Register structures in Julia
2. (optional) Write Patter.map_init() method to sample initial States efficiently
3. Write simulation and density methods:
  - Patter.simulate_obs(state, model, ...)
  - Patter.logpdf_obs(state, model, ...)

Custom workflow

Define a custom State
We will track the location of an animal in three-dimensions:

# Julia code 
struct StateCustom <: Patter.State
  # Map value
  # > This is required for all states & is the value on the map at (x, y)
  map_value::Float64
  # Coordinates
  x::Float64
  y::Float64
  z::Float64
end

Tip

This state is now built into patter as StateXYZ.

Custom workflow

(optional) Define a function to convert initial coordinates to initial states
Here, we simply have to add initial depth values
This is necessary for automated sampling of initial states

# Julia code
function Patter.states_init(state_type::Type{StateCustom}, coords::DataFrame)
    # Add z coordinate
    coords.z = coords.map_value .* rand(nrow(coords))
    return coords
end

Custom workflow

Define a corresponding movement model

# Julia code 
struct ModelMoveCustom{T, U, V, W, X} <: Patter.ModelMove
  # The environment (i.e., map)
  # > This defines the regions within which movements are permitted (i.e., in water)
  map::T
  # Distribution for step lengths
  mobility::U
  dbn_length::V
  # Distribution for headings
  dbn_heading::W
  # Distribution for changes in depth
  dbn_z_delta::X
end

Custom workflow

Define a function to simulate new states using the movement model
This is necessary for sim_path_walk() and pf_filter()

# Julia code 
function Patter.simulate_step(state::StateCustom, model::ModelMoveCustom, t::Int64)
  # Simulate a step length
  length  = rand(model.dbn_length)
  # Simulate a heading
  heading = rand(model.dbn_heading)
  # Calculate new x, y, z coordinates
  x       = state.x + length * cos(heading)
  y       = state.y + length * sin(heading)
  z       = state.z + rand(model.dbn_z_delta)
  # Include the map value
  map_value = Patter.extract(model.map, x, y)
  # Define the state
  StateCustom(map_value, x, y, z)
end

Note

Function arguments are boilerplate. You can edit this code as required for other states. In custom Julia methods, use closure to incorporate additional arguments if required.

Custom workflow

Provide a Patter.logpdf_step() method to compute the unnormalised log-probability of an unrestricted move
This is necessary for pf_smoother_two_filter()

# Julia code 
function Patter.logpdf_step(state_from::StateCustom, state_to::StateCustom,
                            model_move::ModelMoveCustom,
                            t::Int64,
                            length::Float64, heading::Float64)
  # Calculate change in depth
  z_delta = state_to.z - state_from.z
  # Evaluate unnormalised log probability
  logpdf(model_move.dbn_length, length) +
    logpdf(model_move.dbn_heading, heading) +
    logpdf(model_move.dbn_z_delta, z_delta)
end

Custom workflow

(optional) Write an R wrapper to conveniently trial different parametrisations from R:

# R code 
model_move_custom <- function(.mobility = 750.0, 
                              .dbn_length = "truncated(Gamma(1.0, 750.0), upper = 750.0)",
                              .dbn_heading = "Uniform(-pi, pi)",
                              .dbn_z_delta = "Normal(0, 3.5)") {
  glue::glue('ModelMoveCustom(env, {.mobility}, {.dbn_length}, {.dbn_heading}, {.dbn_z_delta});')
}

Custom workflow

Now we’ll define a custom observation model
We imagine we have depth observations for a pelagic species

# Julia code 
struct ModelObsDepthNormalTrunc <: Patter.ModelObs
  sensor_id::Int64
  depth_sigma::Float64
end

Note

This is similar to ModelObsDepthNormalTruncSeabed but the likelihood is normally distributed around the observation, not the seabed. For simplicity, we assume the depth of the seabed is known exactly.

Custom workflow

(optional) Write a Patter.map_init() method for more efficient automated sampling of initial states
This example routine ensures that we only sample from locations where the bathymetric depth is at least as deep at the observed depth

# Julia code 
function map_init(map::Rasters.Raster, 
                  timeline::Vector{DateTime}, 
                  model_move::ModelMove, 
                  dataset::DataFrame,
                  model_obs_type::Type{ModelObsDepthNormalTrunc}, 
                  direction::String = "forward")
    # Check dataset
    check_names(dataset, ["timestamp", "obs", "depth_sigma"])
    # Identify the first depth observation 
    t1  = ifelse(direction == "forward", first(timeline), last(timeline))
    pos = findall(dataset.timestamp .== t1)
    if length(pos) == 0
        return map
    elseif length(pos) != 1
        error("Multiple depth observations recorded at the first time stamp.")
    end 
    depth = dataset.obs[pos]
    # Mask map between limits (false/true -> NaN/true)
    # * map must be >= depth
    msk = (map .>= depth)
    msk = classify(msk, false => missingval(map))
    map = mask(map, with = msk)
    return map
end

map_init (generic function with 2 methods)

Custom workflow

(optional) Write a Patter.simulate_obs() method to simulate new observations
This is necessary for sim_observations()

# Julia code 
# > This simulates s_{t} from s_{t - 1}
function Patter.simulate_obs(state::StateCustom, model::ModelObsDepthNormalTrunc, t::Int64)
  dbn   = truncated(Normal(state.z, model.depth_sigma), 
                    0.0, state.map_value)
  rand(dbn)
end

Custom workflow

Write a Patter.logpdf_obs() method to evaluate the log-probability of the observations
This is required for pf_filter()

# Julia code 
function Patter.logpdf_obs(state::State, model::ModelObsDepthNormalTrunc, t::Int64, obs::Float64)
  dbn   = truncated(Normal(state.z, model.depth_sigma),
                    0.0, state.map_value)
  logpdf(dbn, obs)
end

Custom workflow

Simulate an example movement path:

# Back to R!
state      <- "StateCustom"
mobility   <- 750.0
model_move <- model_move_custom()
path <- sim_path_walk(.map        = map,
                      .timeline   = timeline,
                      .state      = state,
                      .model_move = model_move)

Custom workflow

Simulate an example movement path:

Custom workflow

Simulate observations:

# Run simulation 
 obs <-
   sim_observations(.timeline = timeline,
                    .model_obs = list(ModelObsDepthNormalTrunc = 
                                        data.table(sensor_id = 1L, depth_sigma = 5)))
 arc  <- obs$ModelObsDepthNormalTrunc[[1]]
 yobs <- list(ModelObsDepthNormalTrunc = arc)

 # Visualise simulated observations 
 plot(arc$timestamp, arc$obs * -1,
      col = "royalblue", type = "l",
      ylim = c(-max(arc$obs), 0),
      xlab = "Time stamp", ylab = "Depth (m)")

Custom workflow

Run the algorithms:

 args <- list(.timeline    = timeline,
              .state       = state,
              .yobs        = yobs,
              .model_move  = model_move,
              .n_particle  = 1e5L)

 fwd <- do.call(pf_filter, args, quote = TRUE)
 args$.direction <- "backward"
 bwd <- do.call(pf_filter, args, quote = TRUE)
 smo <- pf_smoother_two_filter()

Custom workflow

Compare simulated depths with particle depths:

plot(smo$states$timestamp, smo$states$z * -1, 
     pch = ".", col = "dimgrey", 
     ylim = c(-max(smo$states$z), 0),
     xlab = "Time stamp", ylab = "Depth (m)")
lines(arc$timestamp, arc$obs * -1, col = "royalblue")

Custom workflow

Compare simulated and reconstructed patterns of space use:

pp <- par(mfrow = c(1, 2))
map_dens(.map = map, 
         .coord = path, 
         .discretise = TRUE, 
         .sigma =  bw.diggle)
map_dens(.map = map, 
         .coord = smo$states, 
         .discretise = TRUE, 
         .sigma =  bw.diggle)
par(pp)

Custom workflow

Parallelisation

In practice, we want to analyse the data for many individuals
You have two options for parallelisation:
1. Iterate over individuals in sequence and use multi-threading in Julia
2. Parallelise over time series in R in a socket cluster and use single-threading in Julia
Forking is not supported with JuliaCall

Parallelisation

Here is some pseudocode to set up a socket cluster

# Initialise Julia
julia_connect(.threads = 1L)

# Set up cluster 
cl <- makeCluster(2L)
clusterExport(cl, ...)
clusterEvalQ(cl, {
  library(patter)
  
  # Iteratively attempt julia_connect()
  # > On Windows, Julia make be locked -> error
  try(julia_connect(), ...)
  
  # Set Julia objects on each core
  JuliaCall:::.julia$cmd("using RCall")
  set_map(...)
  set_seed()
  
  NULL
})

# Run analyses in parallel e.g., via cl_lapply()
# * See .chunk argument for more efficient parallelism
patter::cl_lapply(..., .cl = cl, .fun = function(x) {
  patter_workflow(...)
})

Note

See the patter-eval project for an example where socket clusters were used.

Parallelisation

For large-scale analyses, it is a good idea to gauge (and minimise) requirements
These rack up quickly if you are not careful

Tip

The file_*() helpers (e.g., file_size()) are useful wrappers for interacting with the file system.

Development

Continued improvements and optimisation
Backwards sampling algorithms
patter.workflows package:
- Iterative applications over multiple individuals
- Community-driven conversion routines (e.g., as_actel(), as_glatos())
- Additional helper routines (residency, spatial tools, …)
Community-driven PatterModels.jl accessories package
Other model-based inference packages (e.g., Wahoo.jl)

Issues

Please get in touch with queries, issues & suggestions!
Don’t spend weeks trying to get patter working if you hit a problem—reach out!
See https://github.com/edwardlavender/patter/issues

An introduction to patter

Contents

Contents

Introduction

Introduction

Introduction

Introduction

Team

Institutions

Workshop

Questionnaire

Context

Context

Evolution

Evolution

Evolution

Methodology

Methodology

Methodology

Methodology

Methodology

Methodology

Methodology

Methodology

Methodology

Methodology

Outputs

Why should I use patter?

Why should I use patter?

Why should I use patter?

Why should I use patter?

Why should I use patter?

Set up

Operating systems

Installation

Project set up

Project set up: JULIA options

Project set up: JULIA options

Connect to Julia

Connect to Julia

Workflow

Workflow

Workflow

Workflow

Workflow

Overview

Boilerplate setup (1/7)

Boilerplate setup (1/7)

Boilerplate setup (1/7)

Boilerplate setup (1/7)

Set state (2/7)

Set movement model (3/7)

Set observation model(s) (4/7)

Set observation model(s) (4/7)

Set observation model(s) (4/7)

Run particle filter (5/7)

Run particle smoother (6/7)

Analyse outputs (7/7)

Analyse outputs (7/7)

Analyse outputs (7/7)

Analyse outputs (7/7)

Boilerplate setup

Boilerplate setup

Boilerplate setup

An example workflow with real data

Study system

Study system

Study system

Study system

Study system

Study system

Study system

Study system

Study system

Study system

Study system

Study system

Study system

Study system

Study system

An introduction to `patter`

Why should I use `patter`?

Why should I use `patter`?

Why should I use `patter`?

Why should I use `patter`?

Why should I use `patter`?

Project set up: `JULIA` options

Project set up: `JULIA` options

Connect to `Julia`

Connect to `Julia`